# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD 
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
import matplotlib as mpl
from Spectra.SpectrumA import spectrumAD 

mpl.rc('font',family = 'Times New Roman')

m= 950                           ## Mass in N*s2/m or kg
g= 9.81                          ## Acceleration of gravity in m/s2

adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])

## Model's parameters

k1= 0.55
k2= 0.1
sigAct= 1.5
betaf= 0.5
epsSlip= 0
epsBear= 0
rBear= 1

afs= 0.8
afd= 0.8


# s1p= 3*afs
s1p= 6
# e1p= 2*afd
e1p= 5
# s2p= 9*afs
s2p= 8
# e2p= 10*afd
e2p= 25
# s3p= 35*afs
s3p= 20
# e3p= 55*afd
e3p= 50
pinchX= 0.5
pinchY= 0.01
damage1= 0
damage2= 0
beta= 0

s1n= s1p*10
e1n= np.copy(e1p)
s2n= s2p*10
e2n= np.copy(e2p)
s3n= s3p*10
e3n= np.copy(e3p)

width= 600
height= 1053

testnumber= [13,14,15,16,17,18,19,20,24,25,26,27,28,40,41]

##

datacm= {}
# acceltcm01= {}
# acceltcm02= {}
accelcm= {}
timecm= {}
timecmc= {}

data= {}
datad= {}
accelda= {}
# accelda1= {}
# accelro= {}
accelbase= {}
dispda= {}
dispdacm= {}
timearray= {}
energyarray= {}
energyarraycm= {}
maxaccel= []
minaccel= []
maxaccelf= []
matchindexmin= []
matchindexmax= []
matchdispmin= []
matchdispmax= []
timetotal= []
acceltotal= []
dispnet= {}
dispnetcm= {}
disp1 = {}
disp2 = {}
# timestart= []

for i in range(len(testnumber)):
    datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
    accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]
    timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
    timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]
    disp1['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,5])[indexesb[i]:indexesf[i]+1]
    disp2['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,6])[indexesb[i]:indexesf[i]+1]
    
    
for i in range(len(testnumber)):
    dm = []
 
    for j in range(len(disp1['%d' % (i+1)])):
        d1 = disp1['%d' % (i+1)][j]
        d2 = disp2['%d' % (i+1)][j]
        dm.append((d1+d2)/2)
        
    dispdacm['%d' % (i+1)] = np.array(dm)*adjfactors[i]


data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['8']= np.array(pd.read_csv('20_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['9']= np.array(pd.read_csv('24_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['10']= np.array(pd.read_csv('25_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['11']= np.array(pd.read_csv('26_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['12']= np.array(pd.read_csv('27_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['13']= np.array(pd.read_csv('28_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['14']= np.array(pd.read_csv('40_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['15']= np.array(pd.read_csv('41_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))


for i in range(15):
    accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
    timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
    energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m*(-1),dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(accelcm['%d' % (i+1)]))]
    
cc= 0

for i in range(len(testnumber)):
    for j in range(len(timearray['%d' % (i+1)])):
        if i==0:
            timetotal.append(timearray['%d' % (i+1)][j])
        else:
            timetotal.append(cc+timearray['%d' % (i+1)][j])
        acceltotal.append(accelbase['%d' % (i+1)][j])
    cc= timetotal[-1]
    if i!=len(testnumber)-1:
        for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
            else:
                timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
            acceltotal.append(0)
        cc= timetotal[-1]

#### SDOF Analysis

# for i in range(15):
np.savetxt('TH/Time_Total.txt', timetotal)
np.savetxt('TH/Acceleration_Total.txt', acceltotal)

recordsa= 'TH/Acceleration_Total.txt'
recordst= 'TH/Time_Total.txt'

dts= 0.001953
mass= m/1e6

dur= timetotal[-1]

fh= open('FMdata.tcl', 'w')

fh.write(
    'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
    '\nset M %f;'
    '\nset s1p %f;\nset e1p %f;\nset s2p %f;\nset e2p %f;\nset s3p %f;'
    '\nset e3p %f;\nset pinchX %f;\nset pinchY %f;\nset damage1 %f;'
    '\nset damage2 %f;\nset beta %f;' 
    '\nset s1n %f;\nset e1n %f;\nset s2n %f;\nset e2n %f;\nset s3n %f;\nset e3n %f;'
    '\nset width %d;\nset height %d;'
    '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height, k1, k2, sigAct, betaf, epsSlip, epsBear, rBear))
    # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
fh.close()
subprocess.call('OpenSees MDOF_CalibrationMH_Grav.tcl', shell=True)

####

sdofd= {}
sdofa= {}
sdoft= {}
sdoff= {}
sdofe= {}

dtt = 0.001953

indexes_Beg = [0,int(42/dtt),int(85/dtt),int(127/dtt),int(170/dtt),int(212/dtt),int(255/dtt),int(298/dtt),int(342/dtt),int(387/dtt),int(431/dtt),int(475/dtt),int(518/dtt),int(562/dtt),int(606/dtt)]
indexes_End = [int(28/dtt),int(71/dtt),int(113/dtt),int(156/dtt),int(198/dtt),int(241/dtt),int(284/dtt),int(328/dtt),int(372/dtt),int(416/dtt),int(460/dtt),int(503/dtt),int(547/dtt),int(591/dtt),324987]
timei= [0,42,85,127,170,212,255,298,342,387,431,475,518,562,606]

for i in range(len(testnumber)):
    sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/DispT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])*0.001
    sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
    sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[0][indexes_Beg[i]:indexes_End[i]])-timei[i]
    sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
    sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]) for j in range(len(sdoff['%d' % (i+1)]))]



## Plotting

testnum2= [2,3,4,5,6,7,8,9,11,12,13,14,15,17,18]
intensity= [0.17,0.17,0.19,0.20,0.20,0.19,0.20,0.40,0.29,0.29,0.34,0.39,0.46,0.51,0.52]
location1= [2500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500]
location2= [7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500]
location3= [3.5,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2]
subplotnum = ['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p']

nr= 3
nc= 5

fig1, ax1 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        ax1[i,j].plot(timecmc['%d' % count], energyarraycm['%d' % count], color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax1[i,j].plot(sdoft['%d' % count], sdofe['%d' % count], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax1[i,j].set_xlabel('Time (s)')
            # ax1[i,j].set_xticks([0,10,20,30])
            # ax1[i,j].set_xlim([0,30])
        if count==1 or count==6 or count==11:
            ax1[i,j].set_ylabel('Absorbed Energy (J)')
            # ax1[i,j].set_yticks([0,800,1600,2400,3200])
            # ax1[i,j].set_ylim([0,3200])
        ax1[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, location1[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 9,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax1[i,j].set_xticks([0,10,20,30,40])
        ax1[i,j].set_xlim([0,40])
        
        ax1[i,j].annotate('%s)' % subplotnum[count-1], xy=(35,3700),fontsize=14)
        
        if count==1 or count==2 or count==3 or count==4 or count==5 or count==6 or count==7:
            ax1[i,j].set_yticks([0,1000,2000,3000,4000])
            ax1[i,j].set_ylim([0,4000])
            #ax1[i,j].set_yticks([0,100,200,300,400,500,600])
            #ax1[i,j].set_ylim([0,600])
        else:
            ax1[i,j].set_yticks([0,1000,2000,3000,4000])
            ax1[i,j].set_ylim([0,4000])
        
        ax1[i,j].set_facecolor('#DEE2E6')
        ax1[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        if count==1:
            ax1[i,j].legend(loc=2)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig1.tight_layout()        
#fig1.savefig('MDOF_Absorbed_Energy.png', bbox_inches= 'tight',dpi=300)
fig1.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure10.tiff', bbox_inches= 'tight',dpi=300)


fig2, ax2 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        ax2[i,j].plot(dispdacm['%d' % count]*(-1), accelcm['%d' % count]*m, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax2[i,j].plot(sdofd['%d' % count]*1000, sdoff['%d' % count], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax2[i,j].set_xlabel('Displacement (mm)')
        if count==1 or count==6 or count==11:
            ax2[i,j].set_ylabel('Force (N)')
        ax2[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (-98, location2[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 9,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax2[i,j].set_xticks([-100,-75,-50,-25,0,25,50,75,100])
        ax2[i,j].set_yticks([-10000,-7500,-5000,-2500,0,2500,5000,7500,10000])
        ax2[i,j].set_xlim([-100,100])
        ax2[i,j].set_ylim([-10000,10000])
        
        ax2[i,j].annotate('%s)' % subplotnum[count-1], xy=(75,8000),fontsize=14)
        
        ax2[i,j].set_facecolor('#DEE2E6')
        ax2[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        if count==1:
            ax2[i,j].legend(loc=4)
        
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig2.tight_layout()
#fig2.savefig('MDOF_Hysteretic_Response.png', bbox_inches= 'tight',dpi=300)        
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig2.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure11.tiff', bbox_inches= 'tight',dpi=300)

fig3, ax3 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        pcrdE = np.argmax(abs(dispdacm['%d' % count]))
        pcrdN = np.argmax(abs(sdofd['%d' % count]))
        
        ax3[i,j].plot(timecmc['%d' % count], dispdacm['%d' % count], alpha= 1, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax3[i,j].plot(sdoft['%d' % count], sdofd['%d' % count]*1000, alpha= 0.75, color= '#FF8787', lw= 1.5, label= 'Numerical model')
        ax3[i,j].plot(timecmc['%d' % count][pcrdE],dispdacm['%d' % count][pcrdE],'o',mec='limegreen',mfc='none')
        ax3[i,j].plot(sdoft['%d' % count][pcrdN],sdofd['%d' % count][pcrdN]*1000,'o',mec='#FF8787',mfc='none')
        
        
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax3[i,j].set_xlabel('Time (s)')
        if count==1 or count==6 or count==11:
            ax3[i,j].set_ylabel('Relative Displacement (mm)')
        ax3[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, 90), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 9,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax3[i,j].set_xticks([0,10,20,30])
        ax3[i,j].set_yticks([-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120])
        ax3[i,j].set_xlim([0,30])
        ax3[i,j].set_ylim([-120,120])
        ax3[i,j].set_facecolor('#DEE2E6')
        ax3[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        if count==1:
            ax3[i,j].legend(loc=4)
        ax3[i,j].annotate('%s)' % subplotnum[count-1], xy=(26,97),fontsize=14)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig3.tight_layout()    
#fig3.savefig('MDOF_Relative_Displacement.png', bbox_inches= 'tight',dpi=300)    
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig3.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure12.tiff', bbox_inches= 'tight',dpi=300)

dispmaxE = np.zeros(len(testnumber))
dispmaxN = np.zeros(len(testnumber))

VmaxE = np.zeros(len(testnumber))
VmaxN = np.zeros(len(testnumber))

ratioS2 = [0.6406254468118061,
0.5631207684024964,
0.5982978487537575,
0.5784487347993631,
0.5448047662530605,
0.5272388510824084,
0.5111520596351975,
1.010402694941875,
0.9865847178567868,
0.8809160044546751,
1.047752820467377,
1.0715236408192121,
1.265735924522484,
1.0428579444215016,
1.0407122366162092]

ratioS2V = [1.3152734236888626,
1.252698109663533,
1.3460826893991475,
1.2748593211011638,
1.3007989702302445,
1.2942772102594253,
1.2675454074233086,
1.2909191701902902,
1.1729457235508713,
1.4201618385151407,
1.4484353604472766,
1.368423153566245,
1.6784682056730738,
1.2919167868132813,
1.5419549678869011]

np.mean(ratioS2V)

testnum3 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]

len(ratioS2)

for i in range(len(testnumber)):
    dispmaxE[i] = max(abs(dispdacm['%d' % (i+1)]))
    dispmaxN[i] = max(abs(sdofd['%d' % (i+1)]))*1000
    VmaxE[i] = max(abs(accelcm['%d' % (i+1)]*m))
    VmaxN[i] = max(abs(sdoff['%d' % (i+1)]))
    



fig6, ax6 = plt.subplots(1,2, figsize= (20,8))
ax6[0].plot(testnum3,dispmaxE/dispmaxN,'o',color='#2D708EFF',label='Specimen 1')
ax6[0].plot(testnum3,ratioS2,'o',color='#440154FF',label='Specimen 2')
ax6[0].plot((0,16),(1,1),ls='--',color='k')
ax6[0].set_ylim(0.25,1.4)
ax6[0].set_xlim(0,16)
ax6[0].set_xticks(testnum3)
ax6[0].set_facecolor('#DEE2E6')
ax6[0].set_xticklabels(['2','3','4','5','6','7','8','9','11','12','13','14','15','17','18'],fontsize=14)
ax6[0].set_ylabel(r'$PCRD_{exp}/PCRD_{num}$',fontsize=14)
ax6[0].set_xlabel('Test',fontsize=14)
ax6[0].grid(color= 'black', ls= '--', lw= 0.5)
ax6[0].legend(loc=2,fontsize=14)
ax6[0].annotate('a)', xy=(15.2,0.3),fontsize=16)

ax6[1].plot(testnum3,VmaxE/VmaxN,'o',color='#2D708EFF',label='Specimen 1')
ax6[1].plot(testnum3,ratioS2V,'o',color='#440154FF',label='Specimen 2')
ax6[1].plot((0,16),(1,1),ls='--',color='k')
ax6[1].set_ylim(0.8,2.4)
ax6[1].set_xlim(0,16)
ax6[1].set_xticks(testnum3)
ax6[1].set_facecolor('#DEE2E6')
ax6[1].set_xticklabels(['2','3','4','5','6','7','8','9','11','12','13','14','15','17','18'],fontsize=14)
ax6[1].set_ylabel(r'$Vb_{exp}/Vb_{num}$',fontsize=14)
ax6[1].set_xlabel('Test',fontsize=14)
ax6[1].grid(color= 'black', ls= '--', lw= 0.5)
ax6[1].legend(loc=2,fontsize=14)
ax6[1].annotate('b)', xy=(15.2,0.85),fontsize=16)

fig6.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure18.tiff',dpi=300)


'''
Ta = np.linspace(0,1,100)

SaE = {}
SdE = {}
SaN = {}
SdN = {}

for i in range(len(testnumber)):
    SaE['%d' % (i+1)],SdE['%d' % (i+1)] = spectrumAD(timecmc['%d' % (i+1)],accelcm['%d' % (i+1)],Ta)
    SaN['%d' % (i+1)],SdN['%d' % (i+1)] = spectrumAD(sdoft['%d' % (i+1)],sdofa['%d' % (i+1)]/1000,Ta)


fig5, ax5 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr):
    for j in range(nc):
        ns= nc*i+(j+1)
        ax5[i,j].set_facecolor('#DEE2E6')
        ax5[i,j].grid(color= 'black', ls= '--', lw= 0.5)
        #ax4[i,j].set_xticks([0,10,20,30])
        #ax4[i,j].set_yticks([0,100,200,300,400,500,600,700,800])
        #ax4[i,j].set_xlim([0,30])
        #ax4[i,j].set_ylim([0,800])
        ax5[i,j].plot(Ta,SaE['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax5[i,j].plot(Ta, SaN['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        ax5[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[ns-1],intensity[ns-1]),
                     xy = (0.02, location3[ns-1]), 
                     xycoords = 'data',
                     xytext = (+0, +0), 
                     textcoords = 'offset points', 
                     fontsize = 9,
                     # weight= 'bold',
                     color = 'black',
                     )
        ax5[i,j].annotate('%s)' % subplotnum[ns-1], xy=(0.9,5.5),fontsize=14)
        if ns == 1 or ns == 6 or ns==11:
            ax5[i,j].set_ylabel('Spectral Acceleration (g)')
        if ns == 11 or ns == 12 or ns==13 or ns==14 or ns==15:
            ax5[i,j].set_xlabel('Period (s)')
        if ns == 1:
            ax5[i,j].legend(loc=2)
        
        if ns==13:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,6])
        else:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,6])


fig5.tight_layout()   

#plt.savefig('MDOF_InSpectraAcc.tiff',dpi=300)
fig5.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure13.tiff', bbox_inches= 'tight',dpi=300)
'''